Survival Analysis
9/6/2018
Contents
1 Packagees
library(Vennerable)
library(UpSetR)
library(FDb.InfiniumMethylation.hg19)
library(tidyverse)
library(reshape2)
library(limma)
library(GenomicRanges)
library(DESeq2)
library(pROC)
library(ggplot2)
library(beeswarm)
library(plyr)
library(jyluMisc)
library(glmnet)
library(survival)
library(survMisc)
library(tidyverse)
library(org.Hs.eg.db)
library(DESeq2)
library(jyluMisc)
library(reshape2)
library(pROC)
library(ggplot2)
library(beeswarm)
library(plyr)
library(ggalt)
library(ComplexHeatmap)
library(VennDiagram)
library(limma)
library(Vennerable)
library(survival)
library(ranger)
library(ggplot2)
library(dplyr)
library(survminer)
library(entropy)
library(survMisc)
library(glmnet)
library(rms)
library(grid)
library(gridExtra)
library(maxstat)
library(visNetwork)2 Functions
labeler <- function(x){
x[x==0] <- "Res.WT"
x[x==1] <- "Sens.WT"
x[x==2] <- "Mut"
factor(x)
}
diff.cont <- function(df, cnt){
lev <- levels(factor(cnt))
design <- model.matrix(~0+factor(cnt))
colnames(design) <- c(as.character(lev[1]),as.character(lev[2]))
contrast <- makeContrasts(g0 - g1, levels = design)
fit <- lmFit(t(df), design)
fit <- contrasts.fit(fit, contrast)
fit <- eBayes(fit, 0.01)
tT <- topTable(fit, adjust="fdr", sort.by="B", number=10000)
tT
}
diff_count <- function(df1, df2, f){
df1 <- df1[rownames(df1) %in% rownames(df2),]
f_temp <- df2[rownames(df1), f]
tt <- diff.cont(df1, paste0("g",f_temp))
}
top_sel <- function(tt, frac = 0.1){
rownames(tt)[head(order(tt$adj.P.Val), (dim(tt)[1]*frac))]
}
enricher <- function(DEres, name, gmt){
to.ret <- list()
res <- data.matrix(DEres)
res <- data.frame(res)
res$symbol <- rownames(DEres)
res <- res %>% filter(P.Value < 0.05)
res$P.Value <- res$P.Value * sign(res$logFC)
res <- res %>% dplyr::select(P.Value, symbol)
colnames(res)[1] <- "stat"
res <- res %>% arrange(desc(stat))
rownames(res) <- res$symbol
res$symbol <- NULL
rg <- runGSEA(inputTab = res, gmtFile = gmts[[gmt]])
# rg$Name <- gsub(paste0(name, "_"),"", rg$Name)
rg <- list(rg)
names(rg) <- name
to.ret$rg <- rg
to.ret$enBar <- plotEnrichmentBar(resTab = rg, pCut = .2, ifFDR = FALSE, setName = gmt)
to.ret
}
path.sel <- function(res){
colnames(res)[c(4,6)] <- c("p.up", "p.down")
res %>% dplyr::filter(p.up < 0.2 | p.down < 0.2) -> res
res$Name
}
sel.nut <- function(df, spl, nut, TP53, lambda, alpha){
md1 <-
glmnet(
data.matrix(df),
as.factor(TP53),
family = "binomial",
lambda = lambda,
alpha = alpha,
intercept=FALSE
)
cf <- coef(md1)
imp <- colnames(df)[cf@i]
imp.df <- df[TP53 == 1,imp]
md1 <-
glmnet(
data.matrix(imp.df),
as.factor(nut[TP53 == 1]),
family = "binomial",
lambda = lambda,
alpha = 1,
intercept = FALSE
)
cf <- coef(md1)
cf
}
sel.rate <- function(x){
sum(x != 0)
}
mean.rate <- function(x){
(mean(x != 0)) * ifelse(mean(x) > 0, 1, -1)
}
se <- function(x){
(sd(x != 0))
}
km <- function(response, time, endpoint, titlePlot = "KM plot", pval = NULL, stat = "median") {
#function for km plot
survS <- data.frame(time = time,
endpoint = endpoint)
if (length(levels(factor(response))) > 5){
if (stat == "maxstat") {
ms <- maxstat.test(Surv(time, endpoint) ~ response,
data = survS,
smethod = "LogRank",
minprop = 0.2,
maxprop = 0.8,
alpha = NULL)
survS$group <- factor(ifelse(response >= ms$estimate, "high", "low"))
} else if (stat == "median") {
med <- median(response, na.rm = TRUE)
survS$group <- factor(ifelse(response < med, "low", "high"))
} else if (stat == "binary") {
survS$group <- factor(response)
}
}else{
survS$group <- factor(response)
}
if (is.null(pval)) pval = TRUE
p <- ggsurvplot(survfit(Surv(time, endpoint) ~ group, data = survS),
data = survS, pval = TRUE, conf.int = TRUE,
ylab = "Fraction", xlab = "Time (years)", title = titlePlot,
ggtheme = theme_bw() + theme(plot.title = element_text(hjust =0.5)))$plot
p
}3 Load Preprocessed Data
# setwd("/Volumes/elihei/Internship/projects/TP53ness_ultimate/codes")
load("data/metadata/patmeta_180504.RData")
load("data/Objects/exprMat.RData")
load("data/Objects/CPS1000.RData")
load("data/Objects/ic50.RData")
load("data/Objects/meth.RData")
load("data/patAnnotation/survival_CPS1000.RData")4 Datasets
vennList <- Venn(list(metadata = patMeta$Patient.ID[!is.na(patMeta$TP53)], RNAseq = colnames(exprMat), CPS1000 = rownames(cps.data), methylation = rownames(meth.df)))
Vennerable::plot(vennList, doWeights = F, type = "ellipses", show = list(Faces = FALSE))
upset(fromList(list(metadata = patMeta$Patient.ID[!is.na(patMeta$TP53)], RNAseq = colnames(exprMat), CPS1000 = rownames(cps.data), methylation = rownames(meth.df), IC50 = rownames(ic50.data))), order.by = "freq", number.angles = 30, point.size = 3.5, line.size = 2, sets.bar.color = "orange")
5 Preprocess Data
patMeta <- data.frame(patMeta)
patMeta <- patMeta[!is.na(patMeta$TP53),]
rownames(patMeta) <- patMeta$Patient.ID
colnames(patMeta) <- make.names(colnames(patMeta))
colnames(cps.data) <- make.names(colnames(cps.data))
drug_meta <- merge(patMeta, cps.data, by = 0)
rownames(drug_meta) <- drug_meta$Patient.ID
exprMat <- data.frame(exprMat)6 Define Cut-off for Nutlin Resistency
my_roc <- roc(predictor = drug_meta$Nutlin.3a, response = drug_meta$TP53)
cut.off <- coords(my_roc, "best", ret = "threshold")
cut.off## [1] 0.8670314plot.roc(
predictor = drug_meta$Nutlin.3a,
x = drug_meta$TP53,
print.auc = TRUE,
auc.polygon = F,
grid = c(0.1, 0.2),
grid.col = c("green", "red"),
max.auc.polygon = TRUE,
auc.polygon.col = "blue",
print.thres = TRUE)
outliers <- rownames(drug_meta)[which(drug_meta$TP53 == 0 & drug_meta$Nutlin.3a >= cut.off)]
drug_meta$nut <- ifelse(rownames(drug_meta) %in% outliers, 1, 0)
beeswarm <- beeswarm(
Nutlin.3a ~ TP53,
data = drug_meta,
method = 'swarm',
pwcol = as.numeric(drug_meta$TP53) - as.numeric(drug_meta$nut),
do.plot = F
)[, c(1, 2, 4)]
colnames(beeswarm) <- c("x", "y", "group")
beeswarm$group <- labeler(beeswarm$group)
beeswarm.plot <- ggplot(beeswarm, aes(x, y)) +
xlab("") +
scale_y_continuous(expression("Nutlin.3a Averaged")) +
geom_boxplot(aes(x, y, group = ifelse(group %in% c("Res.WT", "Sens.WT"), "WT", "Mut")), outlier.shape = NA) +
geom_point(aes(colour = factor(group))) +
scale_x_continuous(
breaks = c(1:2),
labels = c("unmutated", "mutated"),
expand = c(0, 0.3)
) + geom_segment(aes(
x = 0.8,
y = cut.off,
xend = 1.2,
yend = cut.off,
fill = "black"
), linetype = "dashed") +
theme_bw() +
labs(x = "TP53", col = "annotation")## Warning: Ignoring unknown aesthetics: fillbeeswarm.plot
# sample.names <- intersect(rownames(drug_meta), colnames(exprMat))
# feature.df <- cbind(drug_meta[sample.names,], t(exprMat[,sample.names]))
# feature.df <- muvis::data_preproc(feature.df)
# feature.df <- feature.df[,-5183]
# mf <- muvis::min_forest(feature.df)7 Differential Analysis with Limma
drug_nut <- diff_count(cps.data, drug_meta[which(drug_meta$TP53 == 0),], "nut")
drug_TP53 <- diff_count(cps.data, patMeta, "TP53")
meth_nut <- diff_count(meth.df, drug_meta[which(drug_meta$TP53 == 0),], "nut")
meth_TP53 <- diff_count(meth.df, patMeta, "TP53")
expr_nut <- diff_count(t(exprMat), drug_meta[which(drug_meta$TP53 == 0),], "nut")
expr_TP53 <- diff_count(t(exprMat), patMeta, "TP53")
top_sel(drug_nut)## [1] "Nutlin.3a" "RO5963" "Cytarabine" "Fludarabine" "Onalespib"
## [6] "Palbociclib"top_sel(drug_TP53)## [1] "Nutlin.3a" "RO5963" "Fludarabine" "Cytarabine" "Selinexor"
## [6] "PF.3758309"top_meth_nut <- top_sel(meth_nut, 0.05)
top_meth_TP53 <- top_sel(meth_TP53, 0.05)
top_expr_nut <- top_sel(expr_nut, 0.05)
top_expr_TP53 <- top_sel(expr_TP53, 0.05)
top_drug_nut <- top_sel(drug_nut, 0.2)
top_drug_TP53 <- top_sel(drug_TP53, 0.2)
vennList <- Venn(list(methylation_nut = top_meth_nut , methylation_TP53 = top_meth_TP53, expression_nut = top_expr_nut, expression_TP53 = top_expr_TP53))
Vennerable::plot(vennList, doWeights = F, type = "ellipses", show = list(Faces = FALSE))
upset(fromList(list(methylation_nut = top_meth_nut , methylation_TP53 = top_meth_TP53, expression_nut = top_expr_nut, expression_TP53 = top_expr_TP53)), order.by = "freq", number.angles = 30, point.size = 3.5, line.size = 2, sets.bar.color = "blue")
8 Drugs
drug.annotation <- read.csv("data/TP53/targetAnnotation_all.csv", sep = ";")
drug.annotation$nameEMBL2016 <- gsub(x = drug.annotation$nameEMBL2016, pattern = " +", replacement = ".")
drug.annotation_nut <- drug.annotation[drug.annotation$nameCPS1000 %in% top_drug_nut | drug.annotation$nameEMBL2016 %in% top_drug_nut,]
targets_nut <- strsplit(as.character(drug.annotation_nut$target), ", ")
targets_nut <- unlist(targets_nut)
drug.annotation_TP53 <- drug.annotation[drug.annotation$nameCPS1000 %in% top_drug_TP53 | drug.annotation$nameEMBL2016 %in% top_drug_TP53,]
targets_TP53 <- strsplit(as.character(drug.annotation_TP53$target), ", ")
targets_TP53 <- unlist(targets_TP53)
vennList <- Venn(list(drug_nut = targets_nut , drug_TP53 = targets_TP53))
Vennerable::plot(vennList, doWeights = F, show = list(Faces = FALSE))
upset(fromList(list(expression_nut = top_expr_nut, expression_TP53 = top_expr_TP53, drug_nut = targets_nut, drug_TP53 = targets_TP53)), order.by = "freq", number.angles = 30, point.size = 3.5, line.size = 2, sets.bar.color = "blue")
# genes <- Reduce(union, list(expression_nut = top_expr_nut, expression_TP53 = top_expr_TP53, top_meth_nut, top_meth_TP53))
DA.features <- Reduce(union, list(expression_nut = top_expr_nut, expression_TP53 = top_expr_TP53, top_meth_nut, top_meth_TP53, top_drug_nut, top_drug_nut))
# sample.names <- intersect(rownames(drug_meta), colnames(exprMat))
# feature.df <- cbind(drug_meta[sample.names,], t(exprMat[genes,sample.names]))
# feature.df <- muvis::data_preproc(feature.df)
# muvis::min_forest(feature.df[,-c(823:dim(feature.df)[2])])9 Gene Set Enrichment Analysis
gmts <- list(HALLMARK=system.file("extdata","h.all.v5.1.symbols.gmt",package="BloodCancerMultiOmics2017"),
ONCOGENIC=system.file("extdata","c6.all.v5.1.symbols.gmt", package="BloodCancerMultiOmics2017"),
KEGG=system.file("extdata","c2.cp.kegg.v5.1.symbols.gmt", package = "BloodCancerMultiOmics2017"))
DEs <- list()
DEs$expression$nut <- expr_nut
DEs$expression$TP53 <- expr_TP53
DEs$methylation$nut <- meth_nut
DEs$methylation$TP53 <- meth_TP53
pathways <- DEs %>% map(function(x) x %>% map(function(y) names(gmts) %>% map(function(z) enricher(y, z, z))))
pathways <- names(DEs) %>% map(function(x) names(DEs[[x]]) %>% map(function(y) names(gmts) %>% map(function(z) enricher(DEs[[x]][[y]], paste(x, y), z))))
pathways[[1]][[1]] %>% map(function(x) plot(x[["enBar"]]))
## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULLpathways[[1]][[2]] %>% map(function(x) plot(x[["enBar"]]))
## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULLpathways[[2]][[1]] %>% map(function(x) plot(x[["enBar"]]))
## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULLpathways[[2]][[2]] %>% map(function(x) plot(x[["enBar"]]))
## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULLresss <- 1:2 %>% map(function(x) 1:2 %>% map(function(y) 1:3 %>% map(function(z) path.sel(pathways[[x]][[y]][[z]]$rg[[1]]))))
ress <- resss %>% map(function(x) x %>% map(function(y) unlist(y)))
vennList <- Venn(list(methylation_nut = ress[[2]][[1]] , methylation_TP53 = ress[[2]][[2]], expression_nut = ress[[1]][[1]], expression_TP53 = ress[[1]][[2]]))
Vennerable::plot(vennList, doWeights = F, type = "ellipses", show = list(Faces = FALSE))
upset(fromList(list(methylation_nut = ress[[2]][[1]] , methylation_TP53 = ress[[2]][[2]], expression_nut = ress[[1]][[1]], expression_TP53 = ress[[1]][[2]])), order.by = "freq", number.angles = 30, point.size = 3.5, line.size = 2, sets.bar.color = "blue")
Reduce(intersect, list(methylation_nut = ress[[2]][[1]] , methylation_TP53 = ress[[2]][[2]], expression_nut = ress[[1]][[1]], expression_TP53 = ress[[1]][[2]]))## [1] "HALLMARK_APOPTOSIS" "HALLMARK_ALLOGRAFT_REJECTION"
## [3] "MTOR_UP.N4.V1_DN" "ESC_V6.5_UP_EARLY.V1_UP"
## [5] "SIRNA_EIF4GI_UP" "KEGG_HEDGEHOG_SIGNALING_PATHWAY"10 Analysis of Elastic Net Coefficients
10.1 With drugs
200 runs
# fdf <- fdf[,!(colnames(fdf) %in% make.names(colnames(cps.data)))]
sample.names <- intersect(rownames(drug_meta), colnames(exprMat))
feature.df <- cbind(drug_meta[sample.names,], t(exprMat[,sample.names]))
name <- rownames(feature.df)
feature.df <- muvis::data_preproc(feature.df)
rownames(feature.df) <- name
# with drugs
# fdf <- feature.df[,-c(5183,5184)]
# lambdas <- (seq(0,2000,10)/100000) * 5 + 0.0001
# cfs <- 1:200 %>% map(function(x) sel.nut(df = fdf[,-c(38, 5177)], nut = feature.df$nut, TP53 = feature.df$TP53, alpha = 0.8, lambda = lambdas[x], spl = 0.75)) %>% map(function(x) data.matrix(x)[,1])
# means <- tapply(unlist(cfs), names(unlist(cfs)), mean)
# ses <- tapply(unlist(cfs), names(unlist(cfs)), sd)
# sel.rates <- tapply(unlist(cfs), names(unlist(cfs)), sel.rate)
# ord <- order(means)
# dat <- data.frame(feature = names(means)[ord], mean = means[names(means)[ord]], sel.rate = sel.rates[names(means)[ord]], ses = ses[names(means)[ord]])
# dat <- dat[dat$mean != 0,]
# #dat <- dat[-grep("ENS",dat$feature),]
# dat %>% arrange(mean) -> dat
# dat$feature <- factor(dat$feature, levels = dat$feature)
# ggplot(dat, aes(x=feature, weight=mean, fill = sel.rate/dim(dat)[1])) + geom_bar() + coord_flip() +labs(y = "mean coefficient", fill = "selection frequency")
# ML.drug.p.features <- as.character(dat$feature)
# write(ML.drug.p.features, "results/ML_drug.txt")
ML.drug.p.features <- read.table("results/ML_drug.txt")
ML.drug.p.features <- as.character(ML.drug.p.features$V1)10.2 Without Drugs
200 runs
# without drugs
# fdf <- fdf[,!(colnames(fdf) %in% make.names(colnames(cps.data)))]
# cfs <- 1:200 %>% map(function(x) sel.nut(df = fdf[,-c(5112)], nut = feature.df$nut, TP53 = feature.df$TP53, alpha = 0.8, lambda = lambdas[x], spl = 0.75)) %>% map(function(x) data.matrix(x)[,1])
# means <- tapply(unlist(cfs), names(unlist(cfs)), mean)
# ses <- tapply(unlist(cfs), names(unlist(cfs)), sd)
# sel.rates <- tapply(unlist(cfs), names(unlist(cfs)), sel.rate)
# ord <- order(means)
# dat <- data.frame(feature = names(means)[ord], mean = means[names(means)[ord]], sel.rate = sel.rates[names(means)[ord]], ses = ses[names(means)[ord]])
# dat <- dat[dat$mean != 0,]
# #dat <- dat[-grep("ENS",dat$feature),]
# dat %>% arrange(mean) -> dat
# dat$feature <- factor(dat$feature, levels = dat$feature)
# ggplot(dat, aes(x=feature, weight=mean, fill = sel.rate/dim(dat)[1])) + geom_bar() + coord_flip() +labs(y = "mean value", fill = "selection frequency")
#
# ML.drug.n.features <- as.character(dat$feature)
# write(ML.drug.n.features, "results/ML.txt")
ML.drug.n.features <- read.table("results/ML.txt")
ML.drug.n.features <- as.character(ML.drug.n.features$V1)vennList <- Venn(list(DA = DA.features , ML_with_drug = ML.drug.p.features, ML_without_drug = ML.drug.n.features))
Vennerable::plot(vennList, doWeights = F, type = "circles", show = list(Faces = FALSE))
upset(fromList(list(DA = DA.features , ML_with_drug = ML.drug.p.features, ML_without_drug = ML.drug.n.features)), order.by = "freq", number.angles = 30, point.size = 3.5, line.size = 2, sets.bar.color = "blue")
intersect(ML.drug.n.features, DA.features)## [1] "ENSG00000242686" "WNK2" "CCDC191" "CCDC50"
## [5] "ENSG00000259344" "ASH2LP3" "ENSG00000261353" "GRAMD4"
## [9] "ENSG00000256482" "ENSG00000264750" "PTP4A1P3" "IGLJ1"
## [13] "ENSG00000270174" "ENSG00000237346" "ENC1" "ENSG00000265380"
## [17] "DMD" "ZNF334" "MYH15" "CA1"
## [21] "FRRS1" "TNFRSF8" "ENSG00000224376" "B3GNT7"
## [25] "NXPH4" "MS4A1" "CRIP3" "RNU1.85P"
## [29] "ENSG00000272403" "NACAD" "CHRNA3" "CD80"
## [33] "TREML2" "FAM57A" "CPXM1" "TNFSF14"
## [37] "ENSG00000180458"intersect(ML.drug.p.features, DA.features)## [1] "ENSG00000242686" "ASH2LP3" "CCDC191" "MYH15"
## [5] "ENC1" "ENSG00000270174" "ENSG00000259344" "DMD"
## [9] "CCDC50" "ENSG00000237346" "ENSG00000273338" "CERS6"
## [13] "GRAMD4" "CEP295NL" "ENSG00000265380" "NXPH4"
## [17] "MS4A1" "LINC02175" "TNFRSF8" "MIR4432HG"
## [21] "NACAD" "CPXM1" "FAM57A" "ENSG00000180458"
## [25] "Fludarabine" "RO5963" "Onalespib" "Cytarabine"intersect(ML.drug.n.features, ML.drug.p.features)## [1] "gain5q" "ENSG00000242686" "CCDC191" "ENSG00000232527"
## [5] "ENSG00000235833" "CCDC50" "ENSG00000205246" "GOLGA8K"
## [9] "ENSG00000259344" "ASH2LP3" "MID2" "ENSG00000217702"
## [13] "GRAMD4" "ENSG00000213222" "SATB1" "BTBD7P1"
## [17] "HAMP" "ENSG00000270174" "ENSG00000237346" "ENC1"
## [21] "ENSG00000227827" "TRIQK" "TAS1R1" "FAM226B"
## [25] "IGHV3.23" "CELSR3" "ENSG00000265380" "CMTM7"
## [29] "ADAD2" "ENSG00000237797" "ENSG00000199899" "ENSG00000249738"
## [33] "LINC01806" "USP44" "DMD" "LGALS9C"
## [37] "MYH15" "DISP3" "IGHV4.55" "FSIP2"
## [41] "TRBV28" "RBMS2" "TNFRSF8" "TMPRSS11E"
## [45] "EEF1DP4" "NXPH4" "IGLV10.54" "OTUD4P1"
## [49] "IL17REL" "ZNF99" "MS4A1" "MAGI1"
## [53] "NR3C2" "FAM66A" "HMGB3P14" "PTPRG.AS1"
## [57] "treatment" "NACAD" "MARK2P11" "LINC01587"
## [61] "KIF13A" "ACVR2B" "KCNC3" "ENSG00000265661"
## [65] "IGHV4.28" "FAM57A" "CPXM1" "ASTL"
## [69] "ENSG00000272823" "ENSG00000180458"features1 <- Reduce(union, list(DA.features, ML.drug.p.features, ML.drug.n.features))11 Survival Analysis on Important features
surv <- cbind(survT, feature.df[survT$patientID,])
fdf <- surv[surv$TP53 == 1,]
to.rem <- which(fdf$OS == 0)
fdf <- fdf[-to.rem,]
fdf <- fdf[!is.na(fdf$OS) & !is.na(fdf$died) & !is.na(fdf$treatedAfter) & !is.na(fdf$TTT),]
OS <- fdf$OS
died <- fdf$died
TTT <- fdf$TTT
treatedAfter <- fdf$treatedAfter
fdf <- fdf[,intersect(colnames(fdf), features1)]
fdf <- fdf[,-c(1,2,dim(fdf)[2])]
cv.fit <- cv.glmnet(data.matrix(fdf), Surv(as.double(TTT), treatedAfter), family="cox", maxit = 1000)
fit <- glmnet(data.matrix(fdf), Surv(as.double(TTT), treatedAfter), family = "cox", maxit = 1000)
Coefficients <- coef(fit, s = cv.fit$lambda.min)
Active.Index <- which(Coefficients != 0)
Active.Coefficients <- Coefficients[Active.Index]
Active.Index %>% map(function(i) km(as.numeric(fdf[,i]), TTT, treatedAfter, sprintf(paste0("Treatment: ", colnames(fdf)[i])), stat = "median")) -> plotList
grid.arrange(grobs= plotList, ncol=2)